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Abstract Ultracold polar molecules in multilayered systems have been experimentally realized very 
, ^ recently. While experiments study these systems almost exclusively through their chemical reactivity, 

^> the outlook for creating and manipulating exotic few- and many-body physics in dipolar systems is 

fascinating. Here we concentrate on few-body states in a multilayered setup. We exploit the geometry 
' of the interlayer potential to calculate the two- and three-body chains with one molecule in each layer. 

I— I The focus is on dipoles that are aligned at some angle with respect to the layer planes by means of an 

external eletric field. The binding energy and the spatial structure of the bound states are studied in 
Q-i several different ways using analytical approaches. The results are compared to stochastic variational 

calculations and very good agreement is found. We conclude that approximations based on harmonic 
^ oscillator potentials are accurate even for tilted dipoles when the geometry of the potential landscape 

^ is taken into account. 

(N 
> 

T-H The prospect of extending the potential of ultracold quantum gases as a quantum simulation device 

into the realm of long-range interacting systems has received a massive boost with the recent success 
in producing cold polar molecules in their rotational and vibrational ground state [Tl|5J[31lll[51[Sl[71[Hlin] • 
Systems of this type hold great promise for exploration of non-trivial few- and many-body dynamics 
PsJ with both fermions and bosons [TUllllj . Three-dimensional samples with dipoles can unfortunately be 

highly unstable to chemical reactions losses [7] . Trapping the dipoles in low-dimensional geometries have 
therefore been proposed to counteract this problem [12,1 3l [Hl[T51ITB] and recent experiments at JILA 
" , I find large geometric effects that confirm a reduced loss in quasi-2D layers [17' and in three-dimensional 

J> optical lattices |18J. The layered geometry with dipolar molecules has attracted much attention and 

exotic many-body phases such as paired states [T9 l [20 l[2n[22|l23l[24l l25] . interlayer coherence akin to 
ferromagnetism "7^, anisotropic superfluidity ^37], density waves [^f2^[5U[l311l321l33j . and non-trivial 
Fermi liquids ^34,,.35,36 . However, the few-body structure is equally intriguing and bound complexes 
have been found with dipoles perpendicular to the plane in bi- [5711551155 1I401HT] and multilayers [IFIHS] , 
for tilted dipoles in single [42] and bilayer setups [43ll44]. and in one-dimensional tubular geometries 
07,48,49,50,51 . 

In this paper we focus on a few-body problem in a multilayered stack of planes when the dipoles 
are not perpendicular to the layer planes. Here we focus on the case of three adjacent planes, but we 
will comment on the implications for more layeres at the end of the presentation. A schematic view of 
the geometry and the dipoles is shown in Fig. [l] The energy and spatial structure of these two- and 
three-body dipolar complexes is our interest and we will use analytical approximations to the potentials 
to construct a harmonic model that is exactly solvable. However, this requires careful consideration of 
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the geometry of the potential landscape of the molecules. To check how well an exact harmonic model 
can describe the system, we performed full stochastic variational calculations in a novel manner that 
takes the deformation of the potential for non-perpendicular dipoles into account. Our findings clearly 
demonstrate that the harmonic approach is very accurate for intermediate and strong dipole strength 
when the potential geometry is included properly. This holds for most tilting angles and only breaks 
down when the dipoles are oriented parallel to the layer. Our finding imply that away from the weakly 
interacting limit, the harmonic approximation works well in layered dipolar systems also in the case 
where the dipoles are tilted. This can be used in exact TV-body harmonic Hamiltonian studies 
153] of both energetics, thermodynamics, and instabilities in larger systems jl31l45l[5^ with more layers 
which is the case for recent experiments with around 30 layers |17j . 




Fig. 1 Schematic picture of the trilayer geometry. The layer distance is d and the coordinates we use are 
indicated in the bottom right-hand corner. The angle between the layer planes and the externally aligned 
dipole moment is 9 and is taken to lie in the ajz-plane. 



2 Model and results 

The dipolar potential between two dipoles with dipole moment D in adjacent layers seperated by a 
distance d is 

x-^ + y-^+iP- i{x cos(g) + dsinje))^ 
^(-' y) = D (,2+,2 + rf2)5/2 ' (1) 

where 9 is the angle between the dipole direction and the layer plane as shown in Fig. [T] Below we 
will use d as the unit for length, and U = mD^/h^d as the dipolar strength with m the mass of the 
molecules. Energies will then be in units of h'^/md'^. The molecules used in recent experiments [17j are 
KRb which has a maximum dipole of 0.566 Debye. With the setup in Ref. |17| . this gives a maximum 
U ~ 1.2 which is in the weak-coupling limit. However, molecules like LiCs or RbCs which are currently 
being cooled have much larger dipole moments [SS] and increasing U by factor of 10 or more should 
be possible with those species. Here we consider the strict two-dimensional limit where we can neglect 
the width of the layers. This holds when the optical lattice potential producing the layers is deep. 
Corrections to this can easily be taken into account by integrating out the transverse (z-direction) 
degree of freedom with the wave function of the lattice (simply a Gaussian when the lattice is deep). 

For 9 = tt/2 the potential is spherically symmetric, but for all other angles there is a distinct 
geometric structure. In Fig. [2] we show a contour plot of the potential for 9 = ^/4 where a clear 
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asymmetric structure in the x-direction is visible. It is this fact that we can exploit to build an accurate 
analytical approximation to the ground-state wave function of the two- and three-body systems. Notice 
that there will be a three-body bound state for any value of U in this setup with one molecule in each 
of the three layers. This can be understood from the fact that the two-body bound state in the case 
of two layers always exists as discussed in Refs. Adding the third layer and molecule will then 

render the binding energy even smaller, i.e. i?3 < i?2. This is a simple generalization of the proof given 
in Ref. [SS] to the layered dipolar case where / dxdyV{x,y) — 0. 
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Fig. 2 Contour plot of the interlayer potential for 9 = n/4 and U = 1. The coordinates x and y are measured 
in units of d. Notice the strong asymmetry of the potential landscape in the a;-direction. The potential has 
reflection symmetry along the j/-direction. 



To further illustrate the interesting geometrical structure of the potential along the x-direction, 
Fig. [3] displays the potential for y = and various values of 9 between zero and 7r/2. For 9 = 0, there 
are two global minima located on each side of a; = 0, whereas for = 7r/2 we have a single minimum 
at x = and also cylindrical symmetry. In the regime < 9 < 9* ^ 0.955 (cos^ 9* = |), the potential 
has two minima, one is located at a: > and is deep, while the other one has a small minimum for 
X < 0. At the intermediate angle 9c {sin^ 9c = |) the potential is maximally asymmetric since here 
the monopole term vanishes |43| . The monopole term can be deduced by expanding the potential in 
a two-dimensional spherical expansion via terms of the form cos (f) where x = r cos and y = r sin (/) 
with r = \l x"^ -V y"^ (the y ^ —y symmetry of the potential excludes sint/) terms). The potential can 
be expanded using the basis function 1, cos </>, and cos 20 only, these correspond to a monopole, dipole, 
and quadrupole term as discussed in detail in Ref. [43] . 

We want to exploit the geometry of the interlayer potential to determine the energy and structure 
of the bound states. This can be done by considering expansions of the Hamiltonian that use the 
minima of the potential (on the x > side) as the starting point. Denote the minimum position by 
(x, y) = (a, 0), where a is defined through the condition ^^^J^'^ — 0. The equation determining a is 

15a(acos(6') d^\VL{9)f - Zaij? -\- d^) - 6cos(6')(a2 + d^){aco%{9) + dsin(6')) = 0, (2) 

from which one can easily see that a scales as d, e.g. a = daa, where oq is the solution of Eq. ^ for 
d=l. 

The potential can now be written in term of the variable w = x — a centered at the minimum 

, ^^w^ +y^ ^-d^ + a'^ + 2'wa-3(wcos(9)+dsin(9) + acos(9))^ 
V{w,y) = D 2 I w2 I 2,0 — ^572 ' 

If we expand this to second order in w and y we have 

V{w,y) =v{d) +a{d)w'^ + P{d)y^. (4) 



Fig. 3 Interlayer potential for y = and for different values of d for (7 = 1. The two special angles are defined 
through the relations cos^ = | and sin^ O,. — ^. 



where 

a{d) = ^ao, (6) 

1 ■^.COC^fell (^''°0-^)(l + °0-3("°('')+°0'^°°(''))^) ^ ao(2ao-6coi,(a)(aoco3(£)) + 5in(8))) 
„ _ ^ ^ "("O + l) 2(a2 + l)2 

"0 - (a;^ + l)5/2 

o(J\ _ o a _ 3 l+ag-5(ao cos(e)+sin(e))^ 

Pl"j — dJrPo , Po — -2 (q2 + 1)7/2 \l) 

Here we have used the scaling property a = dao to extract the dependence on d in each term. This 
will be convenient below. 



2.1 Harmonic three-body chains 

Consider three adjacent layers with interlayer distance d as shown in Fig.[T] The potential for a three- 
particle chain system with one particle in each layer using the harmonic approximation above is 

F(X1,X2,X3, 2/1,2/2,^3) = + ^{X1-X2^ da^f + ^(y, _ y^f 

+ ^(X2 -x^- daof + ^(2/2 - Vz? + - - 2da^f + 0(^2 - 2/3)'. (8) 

Notice the subtraction by 2ao in the term relating the xi and X3 coordinates. This reflects the expected 
geometrical structure of the bound chain. Focus now on the x-depcndent terms. Introducing the relative 

coordinate set xi—x^ — \/2qix, xi — X2 — \J\lix + \J\<l2x, and ^3 — X2 = ^\J\lix + \J\q_x2, we have 
92.) - ^{\^<lix + \^q2x - daof + ^{~\flqix + \^q2x + da^f 



I D^ap Orin "12 _ 17 D'^ag (j„ \2 i U D^ap 2 i QOfso„2 17 ap , /o„ J„ 



ii^(<Zix-y2dao)2 + 3%agL- (9) 



The ground-state energy in this oscillator potential is simply = J^^/Uoq (^-y/f + with one 

term from each of the one-dimensional oscillator coordinates qix and q2x- The 2/-direction will give 
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Fig. 4 Bound state energy of the three-body chain as function of dipolar strength U for 9 — ^. The dot-dashed 
hne (SVM) is the value obt ained from a numerical stochastic variational calculation. The full line (expansion) 
shows the result of Eq. ( 10 1, while the dotted {E[^]) shows the result of using the wave function of Eq. as 
a variational gues s for t lie real dipolar potential of Eq. ([l]). The long-dashed (osc) is the optimized oscillator 
based on Eq. (|14[|, while the short-dashed {E[^osc]) uses me wave function in Eq. (|13[) as a variational guess. 



an identical contribution with /3o instead of ao and the final result for the ground-state energy of the 
three-chain, E, within this harmonic approximation is 



md 
19 



-E 




voU. 



(10) 



Below we will refer to this result for the energy as 'expansion' since it corresponds to a naive expansion 
of the potential around its minimum position to second order. The energies obtained from Eq. ( 10 1 

are shown in Figs. |4]and[5[ For = 9*, we find = '^^^^^ which gives ao ^ 2.66, /Sq ^ 4.11, and 
Vq ~ —1.47. In the perpeiiaicular case where 6 = 7r/2, we have oq = 0, ag = /3o = 6, and vq — —2. 
The corresponding wave function in the 'expansion' method can be written 



^{qix,q2x,qiy,q2y) = iVexp 
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, d 



/2an 



(11) 



where the normalization constant is iV^ 



- y/aof3oU / tt'^ . This can be used to improve the approxi- 



mation by considering instead the expectation value of the full potential from Eq. ([T]) and calculating 



— ^< >^ • The result of this proceduce is also shown on Figs. 



and 



result from Eq. (To|, El^F] is closer to the full numerical result basea on t 
method (SVM) 



5 In comparison to the 
le stochastic variational 

3]T This implies that a harmonic oscillator wave function is a good variational guess. 
The difference between Eq. (10) and El'F] arises mainly because the kinetic energy is overestimated in 
the naive expansion. 

To investigate further the performance of harmonic approximations to the dipolar potential, we 
now consider an optimization on the two-body interlayer potential where the minimum position a is a 
variational parameter. We consider the energy functional 



E[G] = 



< G\H\G > 



<G\G> ' 

where H is the full Hamiltonian with the real potential from Eq. ^ and 



G = exp ( --( — 



' 2^ 



(12) 



(13) 
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is a two-body wavefunction. Fixing U and 9, we can now minimize E[G] to find the optimal set of a, 
/3, and a. This suggests what the optimal oscillator potential should be by replacing the real dipolar 
potential by a harmonic one of the form 

= ( — d — ) + ^^r^ ' 

which will give the energy i?osc — \ ^5 + This optimized harmonic potential for two particles in 

two adjacent layers can now be used as the starting point for the three-body calculations by replacing 
ao, Pq, and oq in Eq. ^ by a, /3, and a (and dropping the constant term vq). The results of this 
calculation is shown on Figs. [4] and [s] where it is labeled 'osc'. Just like above, one can use the 
three-body wave function obtained from this optimized choice of harmonic oscillator to calculate the 
expectation value of the real dipolar potential, the results of which is labeled 'i?[tf'osc]'- 

From Figs. |4] and Figs, [s] we see that our harmonic approximations are close (expect for the naive 
expansion) for t/ > 2 (6* = 7r/2) and U A [9 = 9*). Below these values, the wave function spread 
out too much to be described reliably by a single gaussian. Clearly, the naive expansion approach 



result from Eq. (10 1 performs the worst, although it appears to do somewhat better for tilted dipoles 
(9 < n/2). However, using the wave function from Eq. ( |Tl| ) as a variational guess is a great improvement 
as demonstrated by the Elfp] results. The optimzed oscillator results 'osc' and the corresponding 
variational guess based on this wave function, -E['f'osc]j are both better than the expansion, and they 
differ only slightly. This is to be expected since the optimized oscillator will take as much of the real 
potential into account as possible using the three (essentially free) parameters a, /3, and a. 



The fact that the result of Eq. (10) is rather poor, yet the wave function Eq. (11) is a good 



variational guess, is interesting. Note that the slope of the exact SVM and Eq. ( 10 ) is almost the same. 
This demonstrates that one can get a very accurate approximation for the energy by merely shifting 
the expansion result. This strategy has been pursued for 9 ~ 7r/2 recently to study the properties of 
chains in multiple layers |45[l54j . The current results show that this approximation can be extended to 
the non-perpendicular case also. 

To further investigate the structure we now look at the probability distribution for the three-body 
chain. The (partial) probability distribution for the first molecules is given by 

F{xi,yi) = / \'P{qi q2y)m '^J.'^ ' )dridr2dr3, (15) 

where we use the delta function to eliminate the integration over the coodinates of molecule 1, ri, 
and similarly for the other molecules. In Fig. [6] we show the probability distributions in each layer 
for 9 — 9^ with U — 5 and U — 15 based on the gaussian wave functions used above (expansion 
method) . This wave function very accurately describes the spatial structure of the full solution. The 
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Fig. 6 Probability densities for tiie molecules in each of the layeres when 8 = 0*. Left side has [7 = 5, while 
the right side has U = 15. 



figure clearly demonstrates the geometry of the dipolar potential with maximum in the middle of the 
center layer and displaced maxima in the two outer layers. This trend will continue for more than 
three layers. Notice also how the extend of the probability shrinks as U increases in response to the 
increased localization in the potential minimum. 

While we have only considered 6 — tt/2 and — 9* explicitly, similar results are obtained for a 
large range of angles 6* > 0. Only close to 6* = do we run into problems since the potential has a 
two-minima structure which is not well described with a single gaussian for each molecule. Here one 
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should presumably use at least two gaussians, one for each minimum. This is a subject for future 
investigations. 

3 Discussion 

Wc have studied polar molecules in a geometry consisting of equally spaced two-dimensional planes 
which is predicted to have a rich collection of few-body states due to the long-range and anisotropic 
dipole-dipole interaction between the molecules. In particular, the long-range character reaches across 
the different layers and provides an attractive force that can hold non-trivial bound states of several 
molecules. Here we have considered a simple case where there are three layeres and one molecule in each 
layer. This means that there are only interlayer interactions to worry about. When the dipolc moment 
of the molecules are aligned perpendicular to the layers, a three-body bound chain is always present. 
While this result is expected to hold also when the dipoles are tilted away from perpendicular, no 
systematics has been considered thus far. Here we extend these results and show that for tilted dipoles 
the three-body chain persists and it can be well approximated by various gaussian wave function 
schemes. This conclusion can be easily extended to chains in setups with more than three layers. 

Our results demonstrate that harmonic approximations to the dipolar interaction are very accurate 
in determining the bound state energy and also the spatial structure of the wave function in the 
intermediate and strong dipolar interaction regime also when the dipole moment is tilted with respect 
to the layers. Perhaps somewhat surprisingly, the harmonic approximation seems to work better around 
the critical angle (defined as the angle where two dipoles on a line will have zero interaction) than 
in the perpendicular case. We thus conclude that harmonic approximations are a viable and accurate 
strategy when studying the rich few-body structure of polar molecules in low-dimensional and/or lattice 
geometries. 
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